/*
 * ReactionPath.cpp
 *
 *  Created on: 2 Aug 2011
 *      Author: Allan
 */

#include "ReactionPath.h"

// C++ includes
#include <fstream>
#include <iomanip>
#include <set>

// GeoReact includes
#include "Assemble.h"
#include "KineticsSolver.h"
#include "EquilibriumSolver/EquilibriumSolver.h"
#include "Utilities.h"
#include "Equilibrator.h"

ReactionPath::ReactionPath() :
T(25.0), P(1.0)
{}

ReactionPath::ReactionPath(const vector<string>& allSpecies, const vector<string>& equilibriumSpecies, const vector<string>& kineticSpecies, const vector<string>& initialConditions) :
T(25.0), P(1.0)
{
	Initialiaze(allSpecies, equilibriumSpecies, kineticSpecies, initialConditions);
}

void ReactionPath::Initialiaze(const vector<string>& allSpecies, const vector<string>& equilibriumSpecies, const vector<string>& kineticSpecies, const vector<string>& initialConditions)
{
	multiphase = AssembleMultiphase(allSpecies);
	
	equilibriumReactions = AssembleReactionSystem(equilibriumSpecies);
	
	kineticReactions = AssembleReactionSystem(kineticSpecies);
	
	this->initialConditions = initialConditions;
	
	initialEquilibriumConditions = ExtractEquilibriumConditions(initialConditions, multiphase, equilibriumReactions);
}

void ReactionPath::SetTemperaturePressure(double T, double P)
{
	this->T = T;
	this->P = P;
}

void ReactionPath::Solve(double tInit, double tEnd, double tDelta, const string& fileName)
{
	n = GuessSpeciesComposition(initialConditions, multiphase, T, P);
	
	SolveInitialEquilibrium();
	
	// Print the initial state of the multiphase system
	cout << "Initial State of the Multiphase System: " << endl << endl;
	
	multiphase.PrintState(T, P, n);
	
	KineticsSolver kineticsSolver(multiphase, equilibriumReactions, kineticReactions);
	
	kineticsSolver.SetControlSolution(1.0E-6);
	
	ofstream file(fileName);
	
	file << scientific << setprecision(8);
	
	file << setw(18) << left << "Time";
	file << setw(18) << left << "DeltaTime";
	BOOST_FOREACH(const string& species, multiphase.GetSpecies()) file << setw(18) << left << species;
	file << endl;
	
	double t = tInit, dt = tDelta;
	
	while(t < tEnd)
	{
		cout << "Progress: " << 100.0 * t/tEnd << endl;
		
		file << setw(18) << t;
		file << setw(18) << dt;
		for(int i = 0; i < n.size(); ++i) file << setw(18) << left << n[i];
		file << endl;
		
		kineticsSolver.Iterate(t, tEnd, dt, T, P, n);
	}
}

void ReactionPath::SolveInitialEquilibrium()
{
	// The number of primary species
	const uint Nj = equilibriumReactions.Products().size();
	
	// Check if the number of initial equilibrium conditions equals the number of primary species
	if(initialEquilibriumConditions.size() == Nj)
	{
		Equilibrator equilibrator(multiphase.GetSpecies(), equilibriumReactions.Reactants(), initialEquilibriumConditions);
	
		equilibrator.SetTemperaturePressure(T, P);
	
		equilibrator.Solve(T, P, n);
	}
	else
		cout << "Warning: The number of initial equilibrium conditions is not enough to determine the initial equilibrium state. Proceeding anyway." << endl;
}
